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1.0  Executive  Summary 


During  the  course  of  this  research  effort  computational  parametric  studies  of  the 
microbubble  drag  reduction  phenomena  were  conducted.  The  effects  of  mixture  density 
variation,  free  stream  turbulence  intensity,  free  stream  velocity,  and  surface  roughness  on 
the  microbubble  drag  reduction  were  studied  using  a  single  phase  model  based  on 
Reynolds-averaged  Navier-Stokes  transport  equations.  Additionally,  predictions  of 
Eulerian  multiphase  model  for  microbubble  laden  flow  were  compared  with  Direct 
Numerical  Simulation  from  the  open  literature. 

Good  agreement  was  achieved  between  the  simulations  with  the  single  phase  model  and 
experimental  data  of  Deutsch  et  al.  (2003).  This  good  agreement  was  observed  for  both 
free  stream  velocity  as  well  as  surface  roughness  effect  studies.  Increased  free  stream 
turbulence  intensity  was  observed  to  result  in  lower  drag  reduction,  and  this  effect  was 
stronger  for  higher  density  ratios  of  water  and  injected  gas.  For  the  same  free  stream 
velocity  increase,  the  drag  reduction  was  higher  for  higher  density  ratio.  For  fixed  gas 
injection  rate,  lower  drag  reduction  was  predict  for  higher  free  stream  velocity,  and 
increased  drag  reduction  was  obtained  with  increased  surface  roughness. 

The  drag  reduction  predicted  by  the  multiphase  model  was  substantially  lower  than  that 
predicted  by  the  Direct  Numerical  Simulation  model  of  Ferrante  and  Elghobashi  (2004). 
However,  gas  volume  fraction  and  turbulent  kinetic  energy  profiles  predicted  by  the 
multiphase  model  were  similar  but  not  identical  to  those  predicted  by  the  DNS  of 
Ferrante  and  Elghobashi  (2004). 
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2.0  Introduction 


An  increase  in  range  or  speed  of  the  U.S.  Navy’s  surface  vessels,  submarines,  underwater 
vehicles,  and  weapons  can  be  achieved  by  reducing  the  skin  friction  drag  of  these  objects. 
Microbubble  drag  reduction  (MDR)  is  a  unique  flow  control  technique  that  employs 
injection  of  gas  into  a  liquid  turbulent  boundary  layer  to  form  microbubbles  that  can 
dramatically  reduce  the  skin  friction  drag.  This  technique,  which  is  able  to  provide  drag 
reductions  of  as  much  as  80%,  offers  great  potential  in  Naval  applications. 

McCormick  and  Bhattacharya  (1973)  reported  the  first  microbubble  drag  reduction 
experiments.  During  the  past  decades,  many  research  efforts  have  been  devoted  to 
microbubble  drag  reduction  (Merkle  and  Deutsch  1992).  The  work  conducted  by 
researchers  in  the  former  Soviet  Union  and  in  the  United  States,  primarily  by  the  Applied 
Research  Laboratory  (ARL)  at  The  Pennsylvania  State  University,  provided  the 
benchmark  in  microbubble  drag  reduction  research.  It  has  been  found  that  there  are  many 
factors  that  influence  microbubble  drag  reduction,  including  air  jet  flow  rate,  injection 
process,  free  stream  velocity,  pore  size,  buoyancy,  and  surface  configuration.  As 
evidenced  by  published  papers  in  the  open  literature,  most  of  the  previous  studies  of 
microbubble  drag  reduction  were  conducted  experimentally.  Due  to  the  complexity  of  the 
microbubble  boundary  layers,  theoretical  investigations  have  fallen  behind  the  progress 
made  by  experimental  studies.  It  is  recognized  that  a  better  understanding  of  the 
microbubble  drag  reduction  mechanism  is  critical  to  its  optimal  performance  with 
minimal  use  of  gas  volume  in  practical  applications. 

In  recent  years,  analytical  computational  modeling  of  microbubble  drag  reduction  has 
been  attempted  by  several  researchers  (Madavan  et  al.  1985;  Marie  1987;  Kim  and 
Cleaver  1995;  Meng  and  Uhlman  1998;  Xu  et  al.  2002;  Kunz  et  al.  2003)  to  reveal  the 
mechanism  of  the  microbubble  drag  reduction  phenomenon.  A  study  by  Legner  (1984) 
proposed  a  simple  stress  model  for  gas  bubble  drag  reduction  and  indicated  that  the  drag 
reduction  was  caused  by  a  combination  of  density  reduction  and  turbulence  modification. 
Meng  and  Uhlman  (1998)  suggested  that  bubble  splitting  was  a  plausible  basic 
mechanism  for  reducing  turbulence  in  a  microbubble-laden  turbulent  boundary  layer. 
These  efforts  made  impressive  progress  toward  the  in-depth  understanding  of  the 
mechanism  from  various  angles.  However,  the  available  theoretical  work  is  not  sufficient 
to  answer  all  the  questions  associated  with  microbubble  drag  reduction.  The  relative 
importance  of  postulated  mechanisms  remains  unclear.  One  issue  of  specific  interest  is 
bubble  breakup  during  microbubble  drag  reduction,  which  can  give  rise  to  turbulence 
attenuation.  It  has  also  been  suggested  that  a  simple  density  effect  is  the  dominant  source 
of  drag  reduction.  Very  recent  experiments  at  ARL  have  also  shown  that  significantly 
more  drag  reduction  can  be  obtained  on  rough  surfaces  than  on  smooth  surfaces  with 
microbubble  drag  reduction.  It  is  not  understood  why  this  phenomenon  occurs. 

In  this  project,  the  Applied  Research  Center  (formerly  Hemispheric  Center  for 
Environmental  Technology)  at  Florida  International  University  has  conducted 
computational  parametric  studies  of  the  microbubble  drag  reduction  phenomena.  The 
objective  of  this  research  was  to  assess  the  roles  of  mixture  density  variation,  free  stream 
turbulence  intensity,  free  stream  velocity,  and  surface  roughness  in  microbubble  drag 
reduction.  Additionally,  predictions  of  the  open  literature  Direct  Numerical  Simulation 
for  microbubble  laden  flow  were  compared  with  Eulerian  multiphase  model  results. 
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3.0  Numerical  Model 


3.1  SINGLE  PHASE  MODEL 

During  this  project,  computational  assessment  of  the  role  of  the  mixture  density, 
turbulence  level,  free  stream  velocity,  and  gas  flow  rate  in  the  microbubble  drag 
reduction  was  carried  out.  To  perform  this  assessment,  a  two-dimensional  computational 
fluid  dynamics  (CFD)  model  of  microbubble-laden  flow  over  a  flat  plate  was  developed. 
The  model  consisted  of  Reynolds-averaged  Navier-Stokes  (RANS)  transport  equations 
and  a  standard  k-co  turbulence  model  (Wilcox  1998)  with  a  low  Reynolds  number 
correction.  This  model  was  designed  to  be  applied  throughout  the  boundary  layer  when 
the  near  wall  mesh  is  sufficiently  fine.  The  RANS  transport  equations  are 


(1) 


d_ 
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In  the  equations  (1)  and  (2)  as  well  as  subsequent  equations  p  and  p  are  the  mixture 
density  and  viscosity  respectively,  which  are  defined  in  (12)  and  (13)  below.  The  last 
term  in  equation  (2),  called  Reynolds  stresses,  is  related  to  mean  velocity  gradients  using 
Boussinesq  hypothesis: 
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The  transport  equations  for  the  turbulence  kinetic  energy  k  and  the  specific  dissipation 
rate  0)  are 
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)  is  the  generation  of  turbulence  kinetic  energy 


due  to  mean  velocity  gradients;  G(0  =  a — Gk  is  the  generation  of  t O,  and  Y k  =  p/3  f  .kco 

k  13 

and  Ym  =  pfifpCQ1  are  the  dissipation  of  k  and  co  due  to  turbulence. 
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The  low  Reynolds  number  correction  is  achieved  by  introducing  damp  coefficient  a  into 
the  turbulent  viscosity  equation  as  shown  below  (Wilcox,  1998): 


fk  =  « 


pk 

03 


(6) 


where  a  =  1 


Re, 


a  -  a 


pk 

P03 


go  +  ^er  /R/i 

1 +  Re,  l Rk  j 


(7) 


Rk=  6; 


<  =  4 !  Pi  =  0-072 


The  rest  of  the  coefficients  used  in  the  above  equations  are  defined  as  follows  (Wilcox, 
1998): 
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Wall  boundary  condition  for  the  turbulence  kinetic  energy  is 

^  =  0  (8) 

on 

where  n  is  the  local  coordinate  normal  to  the  wall. 

Wall  boundary  condition  for  the  specific  dissipation  rate  is  discussed  in  detail  below, 
where  the  effect  of  surface  roughness  on  microbubble  drag  reduction  is  addressed. 

Mixture  density  variation  due  to  microbubbles  was  modeled  by  introducing  CCb  gas 
species  and  using  the  species  transport  model.  In  the  species  transport  model,  the  local 
mass  fraction  of  each  species,  Y„  is  predicted  by  solving  a  convection-diffusion  equation 
for  the  zth  species: 

f  (py/)  +  v(pby,)  =  -vy,+5,  (9) 

at 
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where  is  the  diffusion  flux  of  species  i  and  5,  is  the  rate  of  creation  from  any  sources. 
For  turbulent  flows  such  as  those  considered  here,  the  diffusion  flux  has  the  following 
form: 


f 

PDi,m  + 

v 


Sc,) 


VF, 


(10) 


where  Sc,  is  the  turbulent  Schmidt  number,  (ir  is  the  turbulent  viscosity,  and  Di  m  is  the 
diffusion  coefficient  for  species  i  in  the  mixture. 

The  mixture  density  was  computed  using  the  volume-weighted  mixing  law: 


Pm 


1 


(11) 


The  mixture  viscosity  was  computed  using  the  expression  developed  by  Einstein  (in 
Happel  and  Brenner,  1965): 

1  +  2.5*)  (12) 


where  ///  is  the  viscosity  of  the  liquid  and  (/)  is  the  volumetric  concentration  of  gas. 

The  CO2  gas  was  introduced  as  species  mass  source  in  the  first  layer  of  cells  along  the 
porous  section  of  the  flat  plate  (“Wall  2”  on  Figure  1).  Zero-gradient  condition  for  all 
species  was  used  on  the  flat  plate.  This  model  does  not  capture  the  physics  of  the 
microbubbles;  however,  it  allows  one  to  ascertain  whether  a  simple  mixture  density 
variation  effect  is  the  dominant  source  of  microbubble  drag  reduction.  The  model  was 
solved  using  the  FLUENT  6.2  CFD  solver. 

A  712  mm  by  250  mm  computational  domain  (Figure  1)  with  112  x  64  grid  (Figure  2) 
was  used  to  solve  the  model.  The  height  of  the  computational  domain  was  selected  so  that 
it  would  be  at  least  20  turbulent  boundary  layer  thicknesses,  which,  for  the  present  case, 
was  about  8  =  10  mm.  Wall-normal  clustering  of  cells  was  used  to  resolve  the  boundary 
layer,  and  axial  clustering  was  used  near  the  ends  of  flat  plate  sections  (Wall  1,  Wall  2, 
and  Wall3  in  Figure  1).  Dimensions  of  the  computational  domain  and  boundary 
conditions  are  shown  schematically  in  Figure  1.  Constant  velocity  boundary  condition  is 
used  in  the  domain  inlet,  symmetry  boundary  conditions  are  used  in  the  leading  part  of 
the  domain  as  well  as  in  the  far  field  above  the  flat  plate  to  simulate  experimental  water 
tunnel  conditions,  no  slip  boundary  condition  is  applied  at  the  flat  plate,  and  constant 
pressure  boundary  condition  is  applied  at  the  domain  outlet.  Values  of  k  and  co  at  the 
domain  inlet  were  specified  as  1.2x10"  m“/s  and  1.2x10"  s’  ,  respectively.  The  same 
values  were  used  as  initial  guesses  for  k  and  co.  To  compute  gas  inlet  volumetric  flow 
rates,  the  depth  of  the  domain  was  assumed  to  be  102  mm.  The  dimensions  of  the  domain 
corresponded  to  those  of  the  Merkle-Deutsch  (1992)  flat  plate  experimental  configuration 
to  facilitate  comparison  of  the  results. 
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Figure  1.  Schematic  diagram  of  the  computational  domain. 


Figure  2.  The  112  x  64  computational  grid. 


To  validate  that  the  computational  grid  satisfied  y+  <  1  at  the  first  cell  off  the  wall,  the  y+ 
at  the  first  cell  was  plotted  along  the  flat  plate  for  all  computational  cases.  A  sample  plot 
for  the  case  with  Q/lLoA  =  0.01  and  the  free  stream  velocity  of  10.9  m/s  is  shown  in 
Figure  3.  The  y+  values  along  the  flat  plate  for  other  studied  cases  had  similar  values  and 
are  not  presented  here. 
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In  order  to  validate  the  model  further,  simulated  boundary  layer  velocity  profile  for  the 
case  without  gas  injection  was  compared  to  standard  law-of-the-wall  curves.  Boundary 
layer  velocity  profile  in  the  outlet  plain  of  the  computational  domain  is  plotted  in  Figure 
4  in  inner  variables  and  compared  to  the  standard  curves.  It  is  seen  from  the  figure  that 
the  simulated  profile  is  in  good  agreement  with  standard  curves. 

The  ability  of  our  model  to  predict  the  drag  reduction  is  validated  by  comparing  our 
model  predictions  with  the  experimental  data  of  Merkle  and  Deutsch  (1992)  and  with  the 
ensemble  averaged  multifield  two-fluid  modeling  results  of  Kunz  et  al.  (2003).  Figure  5 
compares  computed  drag  reduction  (integrated  drag  coefficient  over  “Wall  3”  flat  plate 
section  normalized  by  its  single-phase  value)  for  a  number  of  nondimensional  gas 
injection  flow  rates  and  several  free  stream  velocities.  The  figure  shows  that  not  only  the 
results  of  the  current  model  are  close  to  those  of  the  more  advanced  model  of  Kunz  et  al. 
but  also  that  the  current  model  correctly  predicts  smaller  drag  reduction  for  lower  free 
stream  velocity  as  observed  in  the  experiments.  The  results  of  our  model  are  very  close  to 
those  of  Kunz  et  al.  model,  however,  both  models  over  predict  experimentally  observed 
drag  reduction  by  as  much  as  50%  in  some  cases. 


Figure  3.  The  y+  in  the  first  cell  off  the  wall  along  the  flat  plate 
(Wall  1  +  Wall  2  +  Wall3)  for  the  case  with  Q/UJ\  =  0.01  and 
IL=10.9  m/s. 
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Figure  4.  Comparison  of  simulated  boundary  layer  velocity 
profile  with  standard  law  of  the  wall  curves. 


Q/U  A 

Figure  5.  Comparison  of  the  computed  integrated  drag 
coefficient  with  the  experimental  data  of  Merkle  and  Deutsch 
(1992)  and  the  numerical  model  of  Kunz  et  al.  (2003). 
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To  facilitate  comparison  of  simulations  with  experimental  study  of  Deutsch  et  al.  (2003) 
on  the  influence  of  surface  roughness  on  microbubble  drag  reduction  a  slightly  modified 
computational  domain  was  used.  This  domain  is  shown  schematically  in  Figure  6.  The 
dimensions  of  this  domain  were  selected  to  closely  approximate  experimental  conditions 
reported  by  Deutsch  et  al.  (2003).  A  127x64  grid  was  used  for  this  domain.  Constant 
velocity  boundary  condition  was  used  in  the  domain  inlet,  symmetry  boundary  condition 
was  used  in  the  far  field  above  the  flat  plate,  no  slip  boundary  condition  is  applied  at  the 
flat  plate  (wall  1,  wall  2,  wall  3,  and  wall  4),  and  constant  pressure  boundary  condition  is 
applied  at  the  domain  outlet.  Values  of  k  and  co  at  the  domain  inlet  were  specified  as 
0.001  nr/s  and  0.1  s'  respectively,  which  corresponds  to  the  turbulent  intensity  of 
0.24%.  The  same  values  were  used  as  initial  guesses  for  k  and  co.  To  compute  gas  inlet 
volumetric  flow  rates,  the  depth  of  the  domain  was  assumed  to  be  146  mm. 


Surface  roughness  was  included  in  the  model  described  above  using  the  equivalent  sand 
grain  roughness  method.  In  this  method  roughness  height  in  inner  variables  is  defined  as 


kl  =  max 


1.0. 


Pksur 

P 


(13) 


Next,  the  value  of  co  at  the  wall  is  specified  as 


co 


w 


p]4 

p 


CO 


(14) 


where  co+  in  the  laminar  sublayer  is  given  by 


c 


co  =  min 


col 


where  col  = 


100 

kt 


K  <  25 
kl  >  25 


(15) 


Figure  6.  Schematic  diagram  of  the  computational  domain  used 
for  surface  roughness  effect  studies. 
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Simulation  results  were  validated  by  comparing  the  results  for  smooth  wall  cases  with 
experimental  data  of  Deutsch  et  al.  (2003).  This  comparison  is  shown  in  Figure  7,  which 
shows  integrated  drag  coefficient  over  “Wall  4”  flat  plate  section  normalized  by  its 
single-phase  value  as  a  function  of  the  volumetric  flow  rate  of  gas  injection.  It  is  seen 
from  the  figure  that  computed  values  are  very  close  to  the  data  fit  of  experimental  values. 
This  good  agreement  validates  our  model  for  further  studies  of  the  influence  of  surface 
roughness. 


Figure  7.  Comparison  of  the  computed  integrated  drag 
coefficient  for  smooth  wall  with  the  experimental  data  of  Deutsch 
et  al.  (2003). 


3.2  MULTIPHASE  MODEL 


A  multiphase  model  was  used  to  compare  predictions  of  standard  CFD  approach  with 
Direct  Numerical  Simulation  approach  recently  applied  to  model  microbubble  drag 
reduction.  The  Eulerian  multiphase  model  was  selected  as  it  produces  most  accurate 
results  among  other  CFD  multiphase  models.  In  this  model,  a  multiphase  flow  is 
described  as  interpenetrating  continua  of  several  phases,  volume  fractions  represent  the 
space  occupied  by  each  phase,  and  the  laws  of  conservation  of  mass  and  momentum  are 
satisfied  by  each  phase  individually.  Thus,  continuity  and  momentum  equations  are 
solved  for  each  phase.  The  Continuity  equation  is: 


1 

Prq 


l  a 


n 


p= i 


-  777, 


J 


(16) 
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where  p  is  the  volume  averaged  density  of  the  qth  phase  in  the  solution  domain. 

The  momentum  conservation  equation  for  a  phase  q  is: 

3  -  a<iVP  +  V'  T,  +  <*qPq9  +  X (KPC, (vP  -  P,)+  ™pqvpq  -  rhqpvj+ 

k/vv+v-k/wv= 

) 


dt 


+  {p  +  p  +  p 

lift.q  ~  1  vm.q  > 


P= 1 

(17) 


=  /  %  r  2  3 

where  t  q  =  aqpq\yvq +VvqT  j+ aq  Aq — pq  V  vqI  is  the  qth  phase  stress-strain 

V  3  J 

tensor;  K  {v  -  v  )  is  an  interaction  force  between  phases;  Fq  is  an  external  body  force; 
Fm  is  a  lift  force;  and  F  is  a  virtual  mass  force. 

The  lift  force  acting  on  a  secondary  phase  p  in  a  primary  phase  q  is  modeled  after  Drew 
and  Lahey  (1993)  as: 


Fm  =  -0.5 pqap (v  -  vp )x  (V  x  vq ) 


(18) 


The  virtual  mass  force  is  encountered  by  accelerating  secondary  phase  particles  (or 
bubbles)  due  to  the  inertia  of  the  primary  phase  mass.  This  force  is  computed  according 
to  Drew  and  Lahey  (1993)  as: 


(  dnvn  dnvn  A 
0-5-tv A,  199  p  p 


dt  dt 


(19) 


where  q —  =  +  (v„ ■  V Xj)  is  the  phase  material  time  derivative. 

dt  dt  q 

To  model  the  interaction  force  between  phases  the  fluid-fluid  exchange  coefficient  is 
defined  as: 


K 


aqapPpf 


pq 


(20) 


PJl 

where  Tp  =  ^  is  the  particle  relaxation  time,  and /is  the  drag  function  modeled  after 
Schiller  and  Naumann  (1935)  as: 


r  _  CD  Re 
J  24 


(21) 


PaVv-Va\di 


where  Re  =  9  p - - —  is  the  relative  Reynolds  number  and  the  drag  coefficient  C/;  is 


Pq 


computed  as: 
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1 24  (l  +  0. 15  Re0'687  )/Re  Re  <  1000 
[0.44  Re  >  1000 


(22) 


To  satisfy  the  conservation  of  energy  an  enthalpy  equation  is  solved  for  each  phase: 

(<*qPq\)+  v  ■  k PA\ ) =  ~aq  :  VS.  -  V  ■  %  +  Z (<k  +  WpAc  ~  ApKp ) <23) 


where  hq  is  the  specific  enthalpy  of  the  c/h  phase,  qq  is  the  heat  flux,  Qpq  is  the  intensity 
of  heat  exchange  between  the  77th  and  qh  phases,  and  hpq  is  the  interphase  enthalpy. 

For  turbulence  modeling  k-e  turbulence  model  for  each  phase  was  employed.  The 
transport  equations  for  this  model  are: 


k  pqK  ) + v  ■  k  Pep  A ) = v  ■ 


qr  q  q  q 

N 


a. 


Mt,q 


+  k  Gk,q-aqpqeq)+ 


±Klq(ciqk,  -  +±K„(u,  -  Uq)- -£*-Va, 

Z=1  Z=1  (JClal  Z=1  aqaq 


(24) 


and 

4(, 

dty 


q'  q  q ' 


q>q  q~q ) 


V- 

(  u 
n  Mtq  V  £ 

+ 

1 A 
- 1 

q  0  q 

V  U£  J 

1  Kl 

C,eaqGk  q  -Cleaqpq£q  + 


C 


3  e 


-uJ-^-Va,  +fx (it  -<?»)■  —Va, 


Z= I 


Z=1 


al(Tl 


1= 1 


(25) 


where  Clq  =  2  ,  Cql  =  2 


'  Vlq  ^ 


1  +  7/, 


_  Tt,lq 


^ F,lq 


To  account  for  interphase  turbulent  momentum  transfer,  the  turbulent  drag  term  is 
computed  as: 


N  /  \  N  -  \  N 
Z K„-  (v,  ~  V, )- =  £ ZT„  (U,  -  U, )-  2 ilrJ, 
1=1  1=1  1=1 


(26) 


where  Ul  and  U  are  phase-weighted  velocities,  and  Fdr .  is  the  drift  velocity  for  phase  1 


defined  as:  vdl  Iq  = 


— ^Vct, 


VlqVq 


J 


This  model  captures  most  of  the  physical  mechanisms  that  are  known  to  be  important  for 
microbubble  drag  reduction  phenomena  such  as  bubble  dynamics  and  turbulence  effects. 
The  model  was  solved  using  the  FLUENT  6.2  CFD  solver. 

A  194  mm  by  290  mm  computational  domain  (Figure  8)  with  120x120  grid  points 
(Figure  9)  was  used  to  solve  the  model.  The  height  of  the  computational  domain  was 
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selected  so  that  it  would  be  at  least  20  turbulent  boundary  layer  thicknesses,  which,  for 
the  present  case,  was  about  8  =  14.2  mm  at  the  exit  of  the  domain.  Wall-normal 
clustering  of  cells  was  used  to  resolve  the  boundary  layer.  Dimensions  of  the 
computational  domain  and  boundary  conditions  are  shown  schematically  in  Figure  8. 
Inlet  velocity  profile  boundary  condition  is  used  in  the  domain  inlet,  symmetry  boundary 
condition  is  used  in  the  far  field  above  the  flat  plate,  no  slip  boundary  condition  is  applied 
at  the  flat  plate,  and  constant  pressure  boundary  condition  is  applied  at  the  domain  outlet. 
Inlet  velocity  profile  was  generated  in  a  separate  domain  so  that  the  turbulent  boundary 
layer  thickness  is  9.7  mm  with  free  stream  velocity  LL=0.83  m/s.  This  velocity  profile  is 
show  in  Figure  10  in  comparison  with  the  experimental  data  of  DeGraaff  and  Eaton 
(2000).  The  simulated  profile  matches  the  experimental  data  reasonably  well  especially  in 
the  boundary  layer  and  in  close  proximity  to  the  flat  plate.  Values  of  k  and  £  at  the 
domain  inlet  were  specified  as  1.2T0"  m  /s  and  1 .44-10"  m“/s  ,  respectively.  The  same 
values  were  used  as  initial  guesses  for  k  and  £.  The  dimensions  of  the  domain  and  the 
inlet  velocity  profile  correspond  to  those  of  Direct  Numerical  Simulation  of  Ferrante  and 
Elghobashi  (2004)  to  facilitate  comparison  of  the  results. 

To  validate  that  the  computational  grid  satisfied  y+  <  1  at  the  first  cell  off  the  wall,  the  y+ 
at  the  first  cell  was  plotted  along  the  flat  plate  for  all  computational  cases.  A  sample  plot 
for  the  case  with  inlet  bubble  concentration  Cv  =  0.01  is  shown  in  Figure  11.  It  is  seen 
from  this  figure  that  the  condition  y+  <  1  is  satisfied  for  the  entire  surface  of  the  flat  plate. 
The  y+  values  along  the  flat  plate  for  other  studied  cases  had  similar  values. 


U(y) 

inlet 

velocity 

profile 


symmetry 


r 

My 

x 

- > 


pressure 

outlet 


0.29  m 


Wall 


I  0.194  m  | 

Figure  8.  The  computational  domain  used  in  the  multiphase  model. 
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Figure  9.  The  120x120  grid  used  in  the  multiphase  model. 


Figure  10.  Inlet  velocity  profile  used  in  the  multiphase  model. 
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Position  along  the  flat  plate,  m 

Figure  11.  The  y+  in  the  first  cell  of  the  wall  along  the  flat  plate  for 
the  Cv  =  0.01  multiphase  model  case. 
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4.0  Turbulence  Intensity  Effect 


The  effect  of  free  stream  turbulence  intensity  was  studied  by  increasing 
turbulence  kinetic  energy  k  from  1x10"  m7s  to  1x10"'  m'/s“  and  than  to  7x10  “  m  /s  . 
The  corresponding  values  of  turbulence  intensity  are  0.24%,  0.75%,  and  2%.  This  study 
was  performed  using  the  single  phase  model  described  in  Section  3.1.  The  results  of  this 
study  are  shown  in  Figure  12  and  Figure  13.  It  is  seen  from  these  figures  that  for  both 
studied  gas  injection  rates  increasing  free  stream  turbulence  intensity  results  in  lower 
drag  reduction  as  indicated  by  increased  value  of  the  CD/CDQ  .  Also  Figure  13  shows  that 
this  effect  is  stronger  for  higher  density  ratios  of  water  and  injected  gas.  For  example,  for 
the  density  ratio  of  DR=559,  increasing  the  turbulence  intensity  from  0.24%  to  2%  results 
in  the  6.3%  decrease  in  the  drag  reduction,  while  the  same  increase  in  the  turbulence 
intensity  for  the  density  ratio  of  DR=125  results  in  only  4.8%  decrease  in  the  drag 
reduction. 


Q/UA=0.01 

Q/UA=0.02 


0.4- 

0.3- 


0.2  4 


o.i  4 — < — i — • — i — > — i — • — i — t — i — ' — i — 1 — i — 1 — i — • — i — > — i — • — i 

0.0  0.2  0.4  0.6  0.8  1.0  1.2  1.4  1.6  1.8  2.0  2.2 

Turbulent  Intensity,  % 

Figure  12.  Free  stream  turbulent  intensity  effect  on  the  drag 
reduction  for  different  gas  injection  flow  rates. 
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Figure  13.  Free  stream  turbulent  intensity  effect  on  the  drag 
reduction  for  different  density  ratios  of  injected  gas  and  water 
(Q/U„A=0.02). 
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5.0  Free  Stream  Velocity  Effect 


The  effect  of  free  stream  velocity  was  studied  by  increasing  the  free  stream 
velocity  from  4.2  m/s  to  7.5  m/s,  than  to  10.9  m/s  and  finally  to  12.4  m/s.  This  study  was 
performed  using  the  single  phase  model  described  in  Section  3.1  with  computational 
domain  shown  in  Figure  1.  Two  nondimensional  gas  injection  rates  Q/U  ,.A  of  0.01  and 
0.02  were  used  in  this  study,  as  well  as  three  density  ratios  of  5,  125,  and  559.  The  results 
of  the  study  are  shown  in  Figure  14  and  Figure  15.  From  Figure  14  it  is  seen  that  for  both 
studied  nondimensional  gas  injection  rates  the  drag  reduction  increases  with  increasing 
free  stream  velocity  as  indicated  by  decreased  value  of  th eCD/CD0.  For  lower 

nondimensional  gas  injection  rate,  steeper  increase  in  the  drag  reduction  for  free  stream 
velocity  increase  from  4.2  m/s  to  7.5  m/s  is  followed  by  a  more  gradual  increase  for  free 
stream  velocity  increase  from  7.5  m/s  to  12.4  m/s.  While  approximately  linear  increase  in 
the  drag  reduction  is  observed  for  higher  nondimensional  gas  injection  rate.  Figure  15 
indicates  that  the  drag  reduction  is  higher  for  higher  density  ratios.  For  example  for  the 
density  ratio  of  DR=559  the  drag  reduction  increased  by  71%  for  free  stream  velocity 
increase  from  4.2  m/s  to  12.  4  m/s,  while  for  DR=5  the  drag  reduction  only  increased  by 
13%  for  the  same  free  stream  velocity  increase. 

A  comparison  of  simulation  results  with  the  experimental  data  of  Deutsch  et  al. 
(2003)  is  shown  in  Figure  16.  The  computational  domain  shown  in  Figure  6  was  used  for 
these  simulations  as  this  domain  closely  mimics  experimental  conditions  of  Deutsch  et  al. 
(2003).  It  is  seen  from  Figure  16  that  the  modeling  results  are  very  close  to  the 
experimental  data.  Same  as  experiments,  modeling  results  for  fixed  gas  injection  rate 
predict  lower  drag  reduction  for  higher  free  stream  velocity.  It  should  be  noted  that 
dimensional  gas  injection  rate  is  used  in  Figure  16,  while  a  non-dimensional  gas  rate  is 
used  in  Figure  14.  Thus,  for  cases  shown  in  Figure  14,  the  dimensional  gas  injection  rate 
had  to  be  adjusted  for  each  free  stream  velocity  to  maintain  fixed  non-dimensional  rate. 
Consequently,  different  behavior  of  drag  reduction  is  observed  depending  on  if 
normalization  is  used  or  not. 
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Figure  16.  Comparison  of  numerical  simulations  with  different 
free  stream  velocities  with  experimental  data  of  Deutsch  et  al. 
(2003). 
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6.0  Surface  Roughness  Effect 


The  effect  of  surface  roughness  was  studied  by  changing  roughness  height  of  the  “Wall 
4”  section  of  the  flat  plate  (see  Figure  6),  which  correspond  to  the  drag  balance  surface  in 
Deutsch  et  al.  (2003)  experiments.  Roughness  heights  ks  of  75  (tm,  150  pm  and  300  pm 
were  used  in  the  simulation  cases  with  free  stream  velocity  of  10.9  m/s  and  gas  injection 
flow  rates  of  0.0  m3/s,  0.002  m3/s,  0.004  m3/s,  and  0.006  m3/s.  This  study  was  performed 
using  the  single  phase  model  described  in  Section  3.1. 

To  quantify  the  amount  of  drag  reduction  achieved  by  introduction  of  microbubbles,  the 
integrated  drag  coefficient  over  “Wall  4”  surface  for  the  case  with  bubbles  was 
nondimensionalized  by  its  value  for  the  case  without  bubbles  obtained  for  the  same 
roughness  height.  The  results  are  presented  in  Figure  17.  The  figure  indicates  that 
increased  drag  reduction  is  obtained  with  increased  surface  roughness  for  a  fixed  gas 
injection  rate.  The  same  trend  was  observed  in  experiments  at  the  Applied  Research 
Laboratory  at  Penn  State  (Deutsch  et  al.  2003). 

Comparison  of  current  numerical  simulations  with  the  experimental  data  of  Deutsch  et  al. 
(2003)  is  shown  in  Figure  18.  This  figure  indicates  that  the  best  agreement  of  numerical 
simulations  with  experiments  is  observed  for  the  smooth  surface.  For  all  rough  surfaces 
the  numerical  simulations  over  predict  the  drag  reduction,  however  this  over  prediction 
reduces  with  increasing  gas  injection  flow  rate. 


Q  ,  m3/s 

gas’ 

Figure  17.  Surface  roughness  effect  on  the  microbubble  drag 
reduction. 
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Figure  18.  Comparison  of  numerical  simulations  with 
experimental  data  of  Deutsch  et  al.  (2003):  a)  Smooth  Surface;  b) 
Fine  Rough  Surface;  c)  Medium  Rough  Surface;  d)  Fully  Rough 
Surface. 
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7.0  Comparison  of  Multiphase  CFD  model  and  DNS 


Predictions  of  multiphase  CFD  model  described  in  the  Section  3.2  were  compared  with 
Direct  Numerical  Simulation  results  of  Ferrante  and  Elghobashi  (2004).  The  boundary 
conditions  of  the  two  approaches  were  matched  to  the  best  possible  extent.  In  both 
simulations  the  free  stream  velocity  was  0.83  m/s.  The  turbulent  boundary  layer  thickness 
on  the  inlet  boundary  of  computational  domains  was  9.7  mm  for  the  both  approaches. 
The  length  of  both  computational  domains  was  194  mm. 

Comparison  of  the  gas  volume  fraction  profiles  at  x=l  82  mm  from  the  domain  inlet 
boundary  is  shown  in  Figure  19.  It  is  seen  from  the  figure  that  the  gas  volume  fraction 
profiles  for  the  smallest  studied  gas  concentration  of  Cv  =  0.001  are  quite  similar. 
However,  at  higher  gas  concentrations  DNS  predicts  sharper  peaks  in  volume  fraction 
profiles,  which  are  further  away  from  the  flat  plate  surface  as  compared  to  the  CFD 
multiphase  model  prediction.  The  profiles  obtained  by  the  multiphase  model  presented 
here  were  obtained  with  a  step  function  gas  volume  fraction  boundary  condition  on  the 
inlet  of  the  computational  domain.  The  boundary  condition  consisted  of  constant  gas 
volume  fraction  inside  the  boundary  layer  and  zero  volume  fraction  above  it. 

An  attempt  was  made  to  obtain  a  better  match  in  the  predicted  volume  fraction  profiles 
by  using  gas  volume  fraction  profile  obtained  by  DNS  model  as  a  boundary  condition  on 
the  inlet  of  the  domain  for  the  multiphase  model.  The  results  are  shown  in  Figure  20.  It  is 
seen  from  this  figure  that  the  resulting  gas  volume  fraction  profile  is  similar  to  the  one 
obtained  with  step  function  profile  boundary  condition,  only  peak  volume  fraction  is 
higher.  A  sharper  peak  predicted  by  DNS  was  not  observed  for  this  modified  boundary 
condition. 

Since  in  DNS  model  the  density  of  the  injected  gas  was  assumed  to  be  zero,  an  attempt 
was  made  to  reduce  the  density  of  the  injected  in  our  multiphase  model.  However,  as 
illustrated  in  Figure  21  the  gas  volume  fraction  profile  for  the  low  density  gas  is  not 
getting  closer  to  the  profile  predicted  by  DNS. 
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Figure  19.  Gas  volume  fraction  profiles  at  x=182  mm. 
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Figure  20.  Comparison  of  gas  volume  fraction  profiles  for 
different  boundary  conditions  of  the  multiphase  model. 
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Figure  21.  Comparison  of  gas  volume  fraction  profiles  for 
different  injected  gas  density  in  the  multiphase  model. 


Comparison  of  the  effect  of  the  microbubbles  presence  on  the  turbulent  kinetic  energy 
(TKE)  is  shown  in  Figure  22.  It  is  seen  from  this  figure  that  the  TKE  predicted  by  the 
multiphase  model  is  not  affected  by  microbubbles  presence,  while  DNS  predicts  decrease 
in  the  TKE  peak  as  well  as  shifting  of  this  peak  away  form  the  flat  plate  surface.  In  is 
worth  noting  that  the  multiphase  model  predicts  TKE  peak  location  at  the  same  distance 
from  the  flat  plate  as  DNS  prediction  for  the  case  with  microbubbles. 

Comparison  of  the  drag  reduction  predicted  by  the  two  models  is  shown  in  Figure  23.  On 
this  figure  the  drag  reduction  is  represented  as  a  ratio  of  the  integrated  drag  coefficients 
over  the  flat  plate  surface  with  and  without  bubbles  in  the  flow.  It  is  seen  from  the  figure 
that  the  drag  reduction  predicted  by  the  multiphase  model  is  substantially  lower  than  that 
predicted  by  the  DNS  model.  It  is  also  seen  that  gas  volume  fraction  profile  used  as  the 
boundary  condition  on  the  domain  inlet  for  the  multiphase  model  has  little  effect  on  the 
predicted  drag  reduction.  This  is  indicated  in  the  figure  by  practically  coinciding  lines  for 
the  step  function  profile  and  the  one  obtained  from  DNS  results. 
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Figure  22.  Microbubble  effects  on  the  Turbulent  Kinetic  Energy. 


Figure  23.  Comparison  of  the  drag  reduction  predicted  by  the 
multiphase  CFD  and  DNS  models. 
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Ferrante  and  Elghobashi  (2004)  reported  gas  concentration  as  the  average  volume 
fraction  of  bubbles  in  the  computational  domain.  However,  in  the  beginning  of  their 
computation  all  bubbles  were  randomly  distributed  inside  the  boundary  layer.  This  leads 
to  bubbles  volume  fraction  in  the  boundary  layer  to  be  3.7  times  higher  than  the  average 
volume  fraction  in  the  domain.  We  have  performed  simulations  for  gas  volume  fractions 
in  the  boundary  layer  which  are  3.7  times  higher  than  the  average  volume  fractions 
reported  in  the  DNS  modeling  (e.g.  we  used  Cv  =  0.0037;  0.037;  0.074  and  instead  of  Cv 
=  0.001;  0.01;  and  0.002  respectively).  The  drag  reduction  predicted  in  these  simulations 
is  compared  with  DNS  model  in  Figure  24.  Compared  to  Figure  23,  the  drag  reduction 
obtained  with  higher  volume  gas  fractions  is  closer  to  that  predicted  by  DNS,  but  the 
difference  is  still  as  high  as  7%.  Additionally,  when  these  higher  concentrations  are  used, 
the  gas  volume  fraction  profiles  differ  significantly  from  those  predicted  by  DNS,  as 
shown  in  Figure  25. 


Figure  24.  Comparison  of  the  drag  reduction  predicted  by  the 
multiphase  CFD  model  with  gas  volume  fractions  increased  3.7 
times  and  DNS  model. 
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Figure  25.  Comparison  of  gas  volume  fraction  profiles  predicted 
by  the  multiphase  CFD  model  with  gas  volume  fractions 
increased  3.7  times  and  DNS  model. 
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8.0  Conclusions 


Major  conclusions  from  the  performed  parametric  numerical  simulation  study  are 
summarized  below: 

•  Increased  free  stream  turbulence  intensity  results  in  lower  drag  reduction,  and  this 
effect  is  stronger  for  higher  density  ratios  of  water  and  injected  gas 

•  For  the  same  free  stream  velocity  increase,  the  drag  reduction  is  higher  for  higher 
density  ratio 

•  For  fixed  gas  injection  rate,  lower  drag  reduction  is  predict  for  higher  free  stream 
velocity,  which  is  in  a  good  agreement  with  experimental  observations 

•  For  a  fixed  gas  injection  rate,  increased  drag  reduction  is  obtained  with  increased 
surface  roughness,  which  is  also  in  a  good  agreement  with  experimental 
observations 

•  Good  quantitative  agreement  was  observed  between  the  simulations  with  the 
single  phase  model  and  experimental  data  of  Deutsch  et  al.(2003) 

•  Gas  volume  fraction  and  turbulent  kinetic  energy  profiles  predicted  by  the 
multiphase  model  were  similar  but  not  identical  to  those  predicted  by  the  Direct 
Numerical  Simulation  of  Ferrante  and  Elghobashi  (2004) 

•  The  drag  reduction  predicted  by  the  multiphase  model  was  substantially  lower 
than  that  predicted  by  the  DNS  model  of  Ferrante  and  Elghobashi  (2004) 
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Nomenclature 

Symbols 

A  Area,  nr 

Cd  Integrated  drag  coefficient,  dimensionless 


DR 

ks 

P 

Q 

Re 

t 

Ui 

u’i 

Wt 


^ water  Density  ratio,  dimensionless 

P  gas 


roughness  height,  m 
Pressure,  Pa 
Flow  rate,  m3/s 

Reynolds  number,  dimensionless 
Time,  s 

Free-stream  velocity,  m/s 

mean  velocity  components  (i  =  1,  2,  3),  m/s 

fluctuating  velocity  components  (/  =  1,  2,  3),  m/s 


friction  velocity,  m/s 


u+  =  —  velocity  in  inner  variables,  dimensionless 

UT 

y  wall  normal  coordinate,  m 

y+  _  Pu,y  wa|]  normaj  coordinate  in  inner  variables,  dimensionless 

P 

(l  Dynamic  viscosity,  Pa-s 

p  Density,  kg/m3 

Tw  wall  shear  stress,  Pa 
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